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ABSTRACT 


The variable aperture seismic array (VASA) is 
used to investigate the upper and lower mantle of the 
earth by measuring daT/dA and the azimuth of arriving 
body wave phases. 

To study upper mantle structure a form of coherency 
velocity spectral analysis employing a zero lag cross- 
correlation technique which scans seismic channels in 
slowness, time, and azimuth has been applied to the P 
codas of a number of teleseismic events collected by 
VASA during 1970. This method has revealed systematic 
patterns of arriving energy which may be attributable to 
reflections of the P wave off discontinuities in the 
upper mantle in the depth range of 130 km to 170 km, and 
at 650 km. Additionally, these patterns seem to contain 
information concerning the source of the earthquake. 

VASA and the other University of Alberta array, 
Peace, and the three Canadian network seismic stations, 
EDM, FSJ, and MCC have been used to make dT/dA deter- 
minations at distances of 80° - 95°. Earthquakes from 
Japan, Asia, and South America have small travel time 
station anomalies and phase velocities that are in good 
agreement with the Jeffreys-Bullen tables. Events from 
Tonga and Samoa have phase velocities that are up to 15% 


higher. The inclusion of the azimuthal deviations of 
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South Pacific events at VASA and LASA permits the anomaly 
to be interpreted as a heterogeneous region of high 
velocity in the mantle at the core-mantle interface with 
a surface projection lying northeast of the island of 


Hawaii. 
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CHAPTER I 


THE SEISMIC ARRAY 


1. The Use of the Array 


In recent years the advent of the seismic array 
has provided geophysicists with a powerful tool for the 
investigation of the structure of the earth's interior. 
With the single station one is presented with a single 
record in time at any point along which there may be not 
one but several phases arriving simultaneously which 
cannot be separated and identified in time. With the 
array, however, one has a number of time records distri- 
buted in space so that events which arrive together in 
time may be separated by velocity and azimuth. This is 
because the small dimension of the array (of the order 
of 150 km) in comparison with the source-to-receiver path 
length results in the wave trains appearing at the 
individual sensors of the array being very nearly 
coherent. The data redundancy then lends itself to 
various array data processing schemes one of which will 
be discussed in the second chapter. In addition the 
array as a whole observes a "bundle" of rays which have 
traversed the interior of the earth in a narrow cone from 
the source thus permitting the detailed investigation of 


localized heterogeneities along this path,as will be seen 


- | re 
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in the third chapter. 

An extremely useful parameter of seismic observa- 
tions is the apparent slowness dT/dA = p, where dT/dA is 
the travel time differential with respect to distance, A, 
and p is the wave parameter. The large aperture seismic 
array can directly measure this parameter by determining 
the apparent velocity (dx/dT) which is the speed and 
direction of the wavefronts of the body wave phases 
generated by an earthquake as they cross the array. The 
usefulness of the parameter p arises from the fact that 
the Wiechert-Herglotz integral (Bullen, 1963) which 
requires p as a function of A is a direct method of in- 
verting seismic velocity-distance observations into 
velocity-depth profiles, providing the earth is spherically 
symmetric. Knowledge of the azimuth of the incoming wave- 
front is useful for determining the amount of deviation 
the ray suffers along its path and whether this deviation 
is large enough to rule out the assumption of the earth's 
approximate spherical homogeneity upon which the Wiechert- 
Herglotz method depends. An array having three or more 
sensors (of a single polarization) can measure both the 
apparent velocity and the azimuth of a phase independent 
of one another. 

In order to study localized effects within the 


earth it is necessary to employ body waves since their 
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shorter wavelengths will permit greater resolution than 
the use of surface waves or normal modes or the long 
period diffracted P (Alexander and Phinney, 1966) or ScS 
(Mitchell and Helmberger, 1973) phases. The uniqueness 
of the structure as determined by inversion techniques 
will then depend on assumptions about the nature of low 
velocity channels, the completeness of the seismic data, 
and the proper identification of the incoming phases 
from sensor to sensor in the array. At teleseismic 
distances (A > 30°) the apparent slowness, dT/dA, as 
determined by an array is independent of the origin 
time of the earthquake and is only slightly dependent 
on source effects such as radiation patterns, near source 
structure, and source depth and location since one is 
dealing with a very narrow cone of nearly parallel rays. 
For an array comparable to the ones used in this study 
having an aperture of 150 km this solid cone of rays has 
angular dimensions (d@.do) of 5" by 57' at a distance of 
lee 

The effect of structure beneath the receiver is 
a more serious problem (Iyer and Healy, 1972) and can be 
empirically evaluated by examining events from many 
different azimuths and distances. This should have the 
effect of revealing only gross subreceiver structure 


since the large aperture of the array when convolved 
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with crustal and upper mantle variations will average out 
all but such gross structure (e.g. a dipping Moho). 
Hence, the smaller the array aperture the more 
negligible will be the source dT/dA and azimuth effects 
and the more restricted will be the cone of rays into 
the receiver system making the resolution of structure 
in this cone more detailed. However, the effect of 
crustal structure beneath the receiver will become much 
greater in the case of the small array and may obscure 
structure encountered elsewhere along the path unless 
the effects are carefully studied and accounted for in 
the dT/dA results. Large aperture arrays of dimension 
150-200 km seem to offer good resolution for the investi- 
gation of both upper mantle and lower mantle structure. 
The methods used in this work have been used 
previously by Niazi and Anderson (1965), Johnson (1967, 
4969)geChinnery and.Toksog {1967),..Husebye,et.al.(1971), 
Husebye (1969), and Montalbetti (1971). Hales and 
Herrin (1972) have done a review of travel times and 


adT/dA observations. 


2. Recording Instrumentation and Array Data Preparation 
The principal network for the acquisition of tele- 
seismic data in 1970 was the portable variable aperture 


seismic array (VASA). The advantage of this portable 
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array over the fixed array is that VASA can be reposi- 
tioned to advantageous locations with respect to the 
earthquake epicenters. The 1970 array consisted of five 
stations located in central Alberta (figure 1) arranged 
in two equilateral triangles 160 and 30 km to the side. 
Each station consisted of three Willmore Mark II seismo- 
meters, a WWVB receiver, amplifiers, a multiplexer, and 
an analog-to-digital converter system. Figure 2 taken 
from Kanasewich et al. (1974) shows a block diagram of 
the modified tripartite digital recording gain ranging 
system used in 1973. However, it does contain the 
essential elements used in 1970. Figure 3 shows the 
response of the amplifier and a Willmore Mark II seismo- 
meter. This indicates that teleseismic P waves (in the 
frequency range 0.7 to 1.8 Hz) and S waves (about 0.2 Hz) 
are relatively unmodified by the system. The signal from 
each sensor was sampled at a rate of 12.5 times/sec and 
was written onto a magnetic field tape holding about two 
days of data. The field tapes were edited and the events 
passed to a station master tape. 

A special WWVB timer was designed which operated 
by first producing 96 sec of synthetic WWVB time signal 
in binary pulse code format generated internally from 
the input hour and minute corresponding to the first 


minute encountered in the data block to be timed 
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Figure 1. Map of sections of Alberta and British Columbia 
showing the location of the arrays used in this study. 


The solid circles denote the individual seismic stations. 


Figure l. 
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Figure 2. Block diagram of the wide-band tripartite 


amplifier system. 
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Figure 2. 
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Figure 3. -Amplifier gain of the tripartite system and 
the combined response of the amplifier and a Willmore 


Mark II seismometer. 
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and then by cross-correlating the synthetic time signal 
with the actual VB signal. The maximum cross-correlation 
corresponded to the coincidence of the two proper minute 
marks and the fraction of the sampling interval or lag 
for which this maximum occurred was determined by a 

three point parabolic fit to the correlation values. 
Thus, from the number of lags the block start time could 
be determined (for good signal to noise characteristics 
the precision of the time determination was to within 1/4 
Gia d SUL ned MLOtVeLL e Apr fl tetas case;..206ms) --Clearly, 
the accuracy of the phase times may become critical if 
the array is small. 

Once accurate data block start times were availa- 
ble the events were edited such that the data for each 
station had exactly two minutes of noise before the 
Jeffreys-Bullen (J-B) predicted arrival time for that 
event .at that station....In addition, the three traces 
on tape were time shifted so that they were each syn- 
chronous with the WWVB time signal. The time shift was 
facilitated by a simple linear interpolation procedure 
which incorporated a correction for the interrecord gap 
between blocks of data. At this stage also, the 
polarities of the signals were checked and any altera- 
tions made. The modified data was then written onto a 
final master tape with all stations reporting an event 


following one another and with each station having a 
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header card identifying the source (NOAA) and receiver 
(VASA) parameters. This arrangement proved to be 
extremely convenient for subsequent study of the events. 
The second array identified in figure 2 as PEACE, 
is operated by the University of Alberta around the 
Bennett Dam on the Peace river in British Columbia and 
was installed in 1967 to determine the effects of crus- 
tal loading. Each station consists of one down-hole 
vertical seismometer the signal from which is recorded 
as a frequency modulation of pulses on magnetic tape. 
Absolute time is furnished by a WWVB receiver. The 
Peace array and the Canadian network station FSJ at 
Fort St. James, British Columbia are 220 km apart on 
a great circle path to the earthquake epicenters in the 
Tonga region. An additional two station combination was 
provided by the two Canadian seismic network stations 
EDM at Edmonton, Alberta, and MCC at Mica Creek, British 


Columbia. 
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CHAPTER II 


VELOCITY SPECTRAL STUDIES OF UPPER MANTLE STRUCTURE 


DY Cintroduction 


In this chapter the concept of velocity spectral 
analysis which employs the coherency of the signal 
(Covespa) will be applied to the P codas of a number 
of shallow teleseismic events as detected by the varia- 
ble aperture seismic array in central Alberta during 
1970. The algorithm for the Covespa, using a zero lag 
cross-correlation technique, has been originally formu- 
lated by the geophysical exploration industry and has 
been modified to scan seismic array data not only in 
slowness (inverse velocity) and time but also in azimuth 
providing a "three dimensional" analysis of the data 


which will be called a Covespagram. 


Z. Coherency Velocity Spectral Analysis 


The established ideas of velocity spectral analysis 
as applied to the study of seismic waves employing array 
data Hagebeen particularly successful in the investiga- 
tion of apparent velocities (the velocities at which 
wavefronts appear to traverse the seismic array) and the 
discrimination of multiples from primary reflected energy 


(Schneider and Backus (1968), Taner and Koehler (1969), 
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Davies, Kelly, and Filson (1971)). Davies et al. have 
successfully employed the concept of velocity spectral 
analysis, Vespa, in the investigation of P waves and 
of core arrivals closely spaced in time. Their method 
has been used by Doornbos and Husebye (1972) in conjunc- 
tion with the cross-correlation method of determining 
relative arrival times of seismic waves (this will be 
discussed in the next chapter) to study core phases. 

The method used by Davies et al. involves the 
formation of a beam by delay and summation of the 
seismic traces of an array and the determination of 
the power in the beam over a specified time window 
which is stepped (usually in increments of 1 sec) down 
the resultant record. This process is then repeated for 
different values of slowness always keeping the beam at 
a constant azimuth until a two dimensional plot 
(Vespagram) of power in slowness and time is generated. 
In actual practice the array is pointed in a constant 
azimuth toward the source and steered over a wide range 
of wave slownesses so that the power in the beam of 
arriving-waves in time may be examined. 

Doornbos and Husebye (1972) have pointed out, 
however, that proper interpretation of the Vespagram 
requires knowledge not only of the response of the 


array as a function of slowness but also as a function 
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of time. The factors which determine the time response 
are essentially the signal variation across the array, 
the duration of the signal, and the asymmetry of the 
Signal. As they have discovered, it may be difficult 

to properly attribute energy to a phase which follows 
very closely in time to a dominant arrival and is 
separated only in slowness due to the leakage of energy 
into sidelobes as a result of amplitude variations across 
the array. 

The process incorporated in the Covespa technique 
involves a normalized zero lag cross-correlation which 
minimizes the problem of sidelobe leakage and reduces 
the likelihood of misinterpretation of the Covespagram 
by accepting only the high coherency measures. The 
equation is a generalization to a two dimensional array 


of a one dimensional form given by Montalbetti (1971), 


heaters Cr a eae ae (>,8S) 
2 Ul 
celers.t) = amt tb : ) ee eee ; 


Lhe 2 2 
) tet : L re 


where, CC is the coherency, M is the number of channels 
(or sensors) in the array, k is an incremental integer 
on channel i (i # k), T is the length of the time window, 


th 


and fs is the amplitude of the i channel at time t. 
’ 


i 
The computation starts by inserting appropriate delays 
into the traces corresponding to a particular slowness, 


s, and azimuth, $6. Then, for each time along the records 
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the zero lag cross-correlations of all combinations of 
two stations within a specified time window (one second 
for compressional, P, body waves and four seconds for 
transverse, S, body waves) are computed, normalized to 
unity and summed. Thus for an array of five stations 
the summation would involve 10 cross-correlation func- 
tions. Since the coherency, CC, is normalized to unity 
it will give a value of unity at a certain time and 
slowness if the phases and shapes of the signal within 
the window at all sensors are the same. In practice, 
the range of acceptable coherencies (0.5 < cc < 1.0) is 
set high enough to ensure that only similar signals, 
irrespective of their power,are plotted. 

By using synthetic seismograms as input to the 
Covespa program it has been determined that although 
low sidelobes do exist for the variable aperture seismic 
array (VASA), Signal strength variation causes no change 
in the sidelobe pattern. This is shown in figure 4 for 
synthetic data at a slowness of 6.0 sec/deg. The effect 
of using the high coherency values is shown in figure 5. 
6.0 db down corresponds to a coherency of 0.5 and 20 db 
down corresponds to a coherency of 0.1. Figure 6 
demonstrates that a prolonged wave train at all sensors 
merely results in an extended Covespa pattern. In 
addition to the slowness-time responses the effect of 


azimuth on the Covespagram has been determined. 
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Figure 4. Change in coherency pattern as a result of 
amplitude variation. In the first Covespagram all 
sensors have equal amplitudes from a synthetic event 
with a slowness of 6.0 sec/deg and with an azimuth of 
230510" Lil thessecona Covespagram the signal amplitudes 
at EDM and RM2 are halved. The insert shows the 
approximate positioning of the array elements with 
respect to EDM. The coherency values are 0.9, 0.8,..., 
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Figure 5. The effect of choosing only high coherency 


values compared with the decibel scale. 
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Figure 6. The effect of a prolonged wavetrain being 


directed into the array. 
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Whitcomb (1973) has pointed out that using delay and 
summation (a linear process) the 10% amplitude point 
for LASA is +16°. The Covespa process with the VASA 
array has a very sharp azimuthal response since for a 
synthetic event, the 10% coherence points are +4° from 
the actual event direction (figure 7). Therefore, this 
process can be used in conjunction with VASA as an 
azimuth-slowness-time discriminator. This has been 
readily accomplished by redesigning the simple Covespa 
program to compute and store complete Covespagrams for 
a number of azimuths (in increments of say one degree) 
producing a "three dimensional" Covespagram (see Appendix). 
Then, for each coherency maximum in time the optimum 
slowness and azimuth are determined by parabolic curve 
fitting. Also, for each maximum the optimum time is 
determined by the same method. This procedure is then 


repeated for each maximum of the Covespagram in time. 


3. Observations 


The teleseismic events used in this discussion 
are from the South Pacific region and Asia (the azimuthal 
ranges being 230°- 260° and 330°- 10° respectively) and 
are in the epicentral distance range of 84°- 95°. These 
events are among the 100 events recorded by VASA during 


1970 and were also used in the study of the inhomogeneity 
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Figure 7a. The effect of azimuth on the Covespagram as 
the array is steered in 1° increments on either side of 
the event direction (235°) “The sevent ssliowness a 
6.0 sec/deg. The position of the array sensors with 
respect; to EDM is sshownsin) the mncert. 
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beneath the island of Hawaii at the base of the mantle 
(see Chapter III). All the events used in this chapter 
have a magnitude of 5.1 or greater. 

Each of the following sets of figures contains 
three separate plots which have been aligned so that 
they relate in time. The upper plot shows the informa- 
tion from the azimuth-slowness-time extractor which has 
detected the coherency maxima of the P coda and plotted 
the azimuthal deviation (true bearing - geocentric 
azimuth) and the slowness variation on expanded scales 
as functions of time along the record. The solid lines 
joining the symbols on these plots are significant only 
as a visual aid. Directly beneath this plot is the 
Covespagram derived from the optimum azimuth depicting 
100 seconds of P coda information and below this is the 
beam formed from the optimum azimuth and velocity. 
Observe that not every point in the uppermost plot 
coincides in time and slowness with a coherency contour 
on the Covespagram. The azimuthal response for this 
process is very sharp and because the Covespagram is 
generated for only one particular direction, energy 
arriving at an angle which deviates from that used for 
the Covespagram may be either suppressed or completely 
obliterated depending on the extent of the deviation. 

Figure 8 shows event 64 which is from the Tonga 


region. ‘The energy is constant in slowness and azimuth 
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Figure 8. Covespagram of event 64 from the Tonga 
Islands: A = 89°, magnitude = 528, and depth = 33 km, 
The upper plot shows the slowness and the azimuthal 
deviations of the coherency maxima in time, the middle 
plot shows the Covespagram for one azimuth, and the 
third plot is the beam formed for one velocity and 


azimuth. 
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for the first 15-20 seconds. The slowness and azimuth 
then exhibit systematic variations at 30 seconds and 
again at 70 seconds from the initial P wave arrival. 
The coherency contours of the Covespagram and the wave 
train of the beam formed trace suggest that energy 
arrives in groups of about 30 seconds in length repeat- 
ing every 35 to 40 seconds. This trailing energy 
phenomenon has been seen before by Davies, Kelly, and 
Filson (1971). 

The Covespagram of this event has been extended 
out to 350 seconds (figure 9) in order to better dis- 
cern the individual energy groupings and the repetitive 
nature of their arrival. After approximately the third 
group the organization which seems to be inherent in this 
pattern begins to break down as the contours become less 
pronounced in time. The PP arrival should appear at 
about 215 seconds after the P arrival with a slowness of 
8.0 sec/deg. It is possible to actually detect this 
phase here although it is weak, possessing a coherency 
level of 0.7. This weakness may in fact be due to an 
azimuthal shift suffered by this ray as it travelled 
through the upper and middle mantle regions. In addition 
there appear to be two strong coherencies at about 140 


seconds and 148 seconds after P. 
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Figure 9. Extended Covespagram of event 64 with the 
beam showing the energy groupings. The predicted 


arrivals P and PP are indicated. 
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Event 88 (figure 10) is another Tonga event show- 
ing similar behavior to the previous case (event 64) 
with strong variations commencing at about 20 seconds 
from the initial P arrival. As is readily discerned in 
these figures, the slowness and azimuth deviations seem 
to be related. The Covespa groupings appear to be 
hidden further down the record as a result of strong 
directional variations (approximately 6° to 10°). A 
number of Covespagrams were computed for a sequence of 
azimuths and superimposed. This composite plot then 
revealed the familiar energy groupings repeating with 
a period of 40 to 45 seconds resembling very closely 
event 64. 

Figure 11 shows an event from the New Ireland 
region. Here one sees very well defined coherency 
groupings 25 seconds or so in duration, repeating at 
40 second intervals. Both the slowness and azimuth 
appear to be quite constant except for the azimuthal 
excursion at about 45 seconds. 

Event 74 (figure 12) is from the Hindu Kush region 
(depth =.210 km) and is quite suggestive of energy group- 
ings in spite of very strong slowness and azimuth varia- 
tions. This observation taken with the broad coherency 
contours of the Covespagram may be indicative of exten- 


sive scattering of the wave train. This behavior seems 
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Figure 10. Covespagram of event 88 from the Tonga 


Islands: A = 91°, magnitude = 5.3, and depth = 33 km. 
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Figure ll. Covespagram of event 72 from the New 
Ireland region: <A = 95°, magnitude = 5.5, and depth = 


20) kn: 
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Figure 12. Covespagram of event 74 from the Hindu 
Kusheregion ss Aspe, magnitude = 5.2, and depth = 


210 km. 
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to be typical of nearly all events from the Asian 
region. 

In addition to these examples a number of events 
from the Aleutians, Japan, and South America were 
analyzed by the Covespa process, and revealed very 
little trailing energy in the P coda with the exception 
of one event from Japan (magnitude 6.1) which indicated 
a moderate manifestation of this phenomenon comparable 
to event 88 (magnitude 5.2) in the definition of the 


coherency contours. 


4. Interpretation 


The first modern seismic investigations of upper 
mantle discontinuities were done by Gutenberg (1960), 
who noticed that PP, pP', and pPP phases (the primed P 
denotes a ray that has passed once through the interior 
of the earth, including the core, from the source to the 
receiver) have more or less emergent beginnings which 
arrive early and are probably due to reflections from 
the Moho or other discontinuities beneath the earth's 
surface. , Subsequently, Hoffman, Berg, and Cook (1961) 
using information derived from the reflection of energy 
from large quarry blasts claimed to have discovered 
evidence for discontinuities at 190, 520, and 910 kn. 


Anderson and Toksoz (1963) studied surface waves and 
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found that there does in fact exist a low velocity 
channel, as indicated by models CIT 8 and 9 for the 

New Guinea- Pasadena data, bounded on the top by a 
gradual decrease in velocity commencing at a depth of 
about 50 km and on the bottom by a sharper increase in 
velocity at about 150 km, the region within being formed 
by the existence of a small amount of partial melt in 

the upper mantle material (Anderson and Sammis, 1970). 
The investigation of upper mantle structure using body 
wave phases has also been done by Archambeau, Flinn, and 
Lambert (1969) who employed spectral amplitudes and 
travel times to discover that abnormally low velocities 
exist within a region extending from approximately the 
base of the crust to a depth of about 150 km. Experi- 
mental observations of P'P' and its precursors as an 
indication of upper mantle structure were done by Adams 
(1968), by Engdahl and Flinn (1969), and by Bolt, O'Neill, 
and Qamar (1969). Adams has discovered world-wide 
reflections at a depth of 65 - 70 km and regional reflec- 
tions at depths of 110- 140 km (Siberia) and 159-190 km 
(Western Europe). More recent work employing this method 
has been done by Whitcomb and Anderson (1970) and 
Whitcomb (1971, 1973) who have found, among other dis- 
continuities, strong reflections at 50 km, at 130 km and 


at 630 km. 
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In light of these investigations on upper mantle 
discontinuities it seems natural to associate the 
observed repetitive energy groupings seen in the Coves- 
pagrams as being due to reflection, single or multiple, 
from such velocity transition zones. However, the 
energy pattern which is nearly constant within each 
group must be interpreted on the basis of the source 
function, pP, sP, and any reverberations within the 
crust. From the time duration of such groups it may 
be possible to derive some information concerning the 
Size of these sources. Assuming a shallow seismic dis- 
turbance propagating at a velocity of 2 or 3 km/sec and 
that a coherent group is typically 20- 30 seconds in 
duration, the source length should be less than 40-60 
km for magnitudes of 5.1 to 6.0. This estimate is com- 
parable with observations of surface rupture length of 
many Shallow earthquakes (Bonilla, Figure 3.16, 1970). 

It has been established that for the observed 
events each energy group is of 20 to 30 seconds in 
extent and that the period of repetitions is in the 
range 35 to 45 seconds. With the exception of event 74 
(figure 12), these events are all situated at shallow 
depths (h< 33 km, where 33 km is often referred to as 
"normal depth") and it is a straightforward matter to 
investigate the theoretical time differences by modelling 


using the phase pdpP proposed by Davies, Kelly, and 
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Filson (1971). Figure 13 shows the scheme envisaged 
where the phase pdpP is a ray which travels from 

the source down to the discontinuity at depth 

d, is reflected upward to the surface (which is 

a strong reflector), and then is reflected again, 
travelling on to the receiver. The time differences 
between arrivals for such rays and direct P at a distance 
of 89° has been calculated and is summarized in table l, 
which is for surface focus earthquakes. The variation 
in the time differences between surface focus and normal 
depth earthquakes was generally less than 0.01 sec. The 
model used for these results is the Birch-Wiggins lower 
mantle velocity model coupled with the University of 
Alberta, UA4, upper mantle profile. 

In keeping with the average observed periods of 
repetition for the energy groups the discontinuity at 
150 km seems most preferable. However, it must be 
emphasized that this is only an average since a small 
number of events does not permit regionalization of the 
discontinuity depths. It appears, though, from the 
variation of the periods that there is a strong reflec- 
tor of energy in the depth range 130 km to 170 km 
corresponding to the sharp increase in velocity at the 
bottom of the low velocity channel. The burst of energy 
which was observed beginning at about 140 sec after the 


initial P arrival and was visible on both the Covespagram 
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Figure 13. Ray plot showing the proposed phase pdpP 


for d = 150° km ata distance or 89°. 


45 


NOSHSONb -- NOITLBANSI Lb SWW HLIM 


ae 


(U1 |b) 


Wh + 





HIY1¢ 





Pagure 12. 





VWO2AIGHA -— WO TRUMITTA BMM HT! 
, 6 58MM WTI (AT JAP Mie 


ah 


46 


Table 1. The time differences between the phase pdpP 
and the direct P arrivals for various depths of dis- 


continuity, d, are shown for surface focus earthquakes. 


Depth of Discontinuity Time Difference 


(km) (pdpP-P) (sec) 
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and the summed trace (figure 9), may be indicative of 
p6o50pP. Slight variations in dip between successive 
multiple reflectors probably cause the velocity and 
azimuth variations which begin by being relatively con- 
stant and then assume more pronounced deviations. The 
immediate and strong variations in velocity and azimuth 
of the arrivals from event 74 (figure 12) indicate the 
structural complexity of the Asian region in general. 
Questions concerning the location of the reverberations, 
whether beneath the source or receiver, cannot yet be 
satisfactorily answered and may prove to be a combination 
of the two. This may only be resolved by the inclusion 
of a large amount of data from earthquakes of varying 


magnitudes, depths, azimuths, and distances. 


Ss, Conclusion 


The Covespa process used here presents a clear 
picture of the seismic arrivals in azimuth, slowness, 
and time. It is particularly useful in examining tele- 
seismic events which have well developed P coda for the 
purpose of extracting information concerning the earth- 
quake sonres and upper mantle structure. This method 
seems applicable for most earthquakes of magnitude 
greater than 5.1 whereas the P'P' precursor method 
requires earthquakes of magnitude 6.0 or greater. 


Therefore, far better spatial distribution of events 
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in relation to regions of tectonic activity may now 
prove to be available for the study of upper mantle 
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CHAPTER III 


SEISMIC ARRAY EVIDENCE FOR A LOWER MANTLE HETEROGENEITY 


BENEATH THE ISLAND OF HAWAII 


pap, het a a fon layenb ayy 


In this chapter the detection of the first arrivals 
of body waves by the seismic array (VASA) and the deter- 
mination of the wave parameter p = dT/dA will be discussed. 
Two of the methods used to calculate the wave slowness 
involve the inversion of the relative arrival times of 
the body phases at the individual sensors of the array, 
and the third method involving the application of 
coherency techniques to the P coda of array data has 
essentially been discussed in the preceding chapter. 

The results of the dT/dA measurements will then 
be discussed in the context of the earth's mantle 
velocity structure in general, and the structure near 
the base of the mantle underlying the geologically and 
geophysically interesting Hawaiian islands (Chapter IV) 
in. particular.,. The ability of an array»possessing three 
or more sensors to independently determine both the 
apparent velocity (inverse slowness) and the true 
azimuth of an arriving phase will be used in a detailed 
study of the behavior of body waves passing through this 


region. 
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2. Array Determined dT/dA Methods and dT/dA Observational 


Results 


In 1970 the VASA array in central Alberta recorded 
approximately 100 events, eight being from the Tonga and 
Samoa islands region. Additionally the Peace array and 
stations FSJ, MCC, and EDM recorded a total of 32 Tonga 
earthquakes from 1969 to 1971 having good signal-to-noise 
characteristics and magnitudes greater than 5.3. This 
comprised the raw data set. 

aT/dA calculations from the VASA digital data were 
carried out using three methods. The first two methods 
both involved first arrival phase time determinations 
by a semiautomatic technique employing the cross-correla- 
tion of approximately the first cycle and a half of the 
foetal thee Veleror tall Station pairs.’ Sance all 
station pairs produced redundant time differences this 
permitted checks to be made on the data scatter. For 
the critical»south Pacific events the length of the data 
string to be cross-correlated and its position along the 
record were varied somewhat and the correlation procedure 
was repeated thus providing a check on the stability of 
the time determinations. Time determinations were also 
ilade by hand for all events. 

The first method of inverting the time differences 


to velocities and azimuths was the direct three-station 
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inversion. The array was first broken down into all 
combinations of three stations and the following 
analysis was performed on each triangle. 

Let the station triangle be as shown in figure 14. 


The velocity along leg 1,2 is 


ts (3=1) 


12 = 412/712 


where Ai 15> the Gistance from station | to station 2 and 


Tho is the time difference between the wavefront arrivals 


at the two stations. Similarly for leg 1,3 one has 


C (Se2)). 


13 ~ 4337713 


Then the apparent velocity of the wavefront will be given 


by 

C= Ci sin vz (3-3) 
and also 

ae C13 sin(A + a) (3-4) 


where A is the angle formed by leg 1,2 and the wavefront 
and a is given by the known triangle geometry (angle 


213). Solving (3-3) and (3-4) for C one has 


tan A sin A/cos A 


sin a/ (Cj 4/C}3 - cos a) (3-5) 


giving 
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Figure 14. The three station velocity and azimuth 
determination. The array geometry is given by angle 
a and distances Al» and Aj3: The wavefront enters 
the array at station 1 with velocity C, making an 


angle, A, with the side (12). C andsG are the 


12 1S) 
projections of C along legs (12) and (13). 
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Co COS(90-A) = SIN(A) 
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a COS(A+ «-90) = SIN(A) 


Figure 14. 
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AtALaAme [Sin a/ (Cy 5/C14 - cos a)! (S2egt .< 


Hence C may be found using equation (3-3). 

Thus the azimuth and the apparent velocity may 
be found for station 1. This is then repeated for all 
such triangles. The individual velocities and azimuths 
are weighted according to the relative area of their 
respective triangles. The weighted velocities are then 
averaged and the weighted azimuths are transformed to a 
common point (e.g. the center of the array) and then 
averaged. 

The second technique involves a least square 
method for direct measurement of daT/dA and is described 
by Husebye (1969). A brief outline of the method will 
be presented here. 


Let 

eee eee TUB ae oT Ne es Se (3-7) 
represent the response surface of the relative arrival 
time over an arbitrary array. The {B, } are coefficients 


to be determined by a least squares fitting of the sur- 


face to the data and the {25} are defined as follows: 


(3-8) . 


The {2, } are the transformed coordinates of the array 


station locations ¢ (latitude) and \ (longitude) and the 
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exponents Al and A2 are to be determined by a first order 
expansion of the fitted surface T about the first guesses 
ofent, BAZ. 

Since T is the fitted surface determined from the 
guesses (initially Al=1.0, A2=1.0) expanding about 


Al,A2 to first order gives 


oT e) 


Z = oT : aS = 
Top= Pt (Al a0) yAL + (A2 rare tp) TAD (3-9) 
re fe 
= ss oo a 
Top = it (Al ad 36 9 In¢<-+ (A2—-1.0) yy A ind 
(3-10) 
where Tob are the expected values of the arrival times. 


Equation (3-10) is now treated as a new surface given by 


T=B.2z,+B.az,+...+ Bz. + (Al-1.0)z + (A2-1.0)z 


heen oe ge kk are k+2 
(3= 11) 
where 
Zh 4 = = do ln o (3-12) 
and 
Zz = a xX Lik s— les ake 


kz 
The approximate values 93T/3¢ and 3T/3A are determined 
by differentiating the fitted response surface T given 


by (3-7). Now, a new surface is formed in the manner of 


(3-7): 
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T= BAZ .teBSZ He. SeeeB.A zed B 


idk Udit tae Kee Bye apy tee eS Aya re (8-44) 


where Bay = Al-1.0 and Bia = A2-1.0. This new response 


surface is fitted in a least squares fashion to the data 


giving values for B and Buy from which are derived 


k+1 


the new values for Al and A2, 


2 


Al = Bal FLO" ; (3-15) 


A2 = Bi 42 ae oy to-209 


which are used in place of the first guesses in a new 
iteration. The iterations are repeated until either 
certain preset error criteria are satisfied or until no 
further improvement occurs. For the VASA array the 


response surface 


rem (s0+ biG a + BK + B,g (419) 


proved satisfactory requiring only a single iteration 
to reach a stable solution. Large arrays possessing 
apertures greater than 3° will require a response surface 


of second order in @¢ and }: 


see: OB, oo 2 Sees r eee ig hee rn Eee 


(3-18) . 
An array similar to Husebye's (1969) with an aperture 


of 10° and fifteen randomly placed sensors was tested 
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using (3-18) and required two iterations to reach a 
stable solution. 

The apparent slowness is determined by finding 
the orthogonal comoonents for each station i, aT/3o 5 
and dT/OA5 5 in the north and east directions respectively 
using the best fitting time response surface, T. The 


resultant slowness P; is given by 


- 2 2 
Pp, =7u; + Vi (3-19) 
“a af « 
h = Saas So 
where u a0; =H and Ve hy e, ; 


e. ande being unit vectors pointing north and east 
respectively. 


The azimuth is found using the equation 
tan ie = v;/u; (3-20) . 


In a similar manner dT /ay* may also be found. These 
quantities need not be calculated at the station loca- 
tions but, in general, can be determined for any set of 
points within the array. 

Errors are determined by assuming that errors in 
the various parameters (u;, EAE ee 0 iG «e») are inde- 
pendent Gaussian variables having zero mean and then 
finding the first order perturbations on the parameters 
(Husebye, 1969). Computationally the problem is set up 


in the form of normal linear equations and solved using 
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matrix operations: from equation (3-7) with k = 4 the 


normal equations will be 


L257 = 22424 + BiuZ,Z5 4 B3%212, + BAZ 724 


tZ5T = LZ5Z1 + BitZ525 + B,%252, + B,i2oZ4 

LZ 7 = UZ3Z) oP B,2232Z, + B22323 + B,%Z32Z4 

LZ ,T = LZ ,Zy + BiUZ4Z5 “9 B,2242, + B,U2 42, 
(3-21) 


where the Zz, are defined by (3-8) and the summation runs 
over all points of observation. (3-21) can be put into 


matrix form: 


C=] Bd (3-22) 


where C is the column vector whose elements are the 
left side of (3-21) and B is the column vector formed 
by the coefficients {B, } and Z is the matrix formed by 
the summation products of the {2,}. The desired co- 


efficients {B, } can be found by solving (3-22): 
Bo = C.2 (om 2 3) 2 


The third technique of velocity-azimuth deter- 
mination involves the nonlinear array processor employing 


a coherency measure requiring only the array geometry and 
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the block start times. This method was discussed in 
Chapter IL. 

aT/dA results for the three methods proved to 
be in close agreement being within 0.01 sec/deg for 
those events for which comparison was made (the South 
Pacific and Asian events). Also, it was observed that 
the hand picked first arrival times yielded velocities 
and azimuths that were very nearly the same as those 
obtained using the cross-correlation techniques although 
with the latter methods the data scatter was considerably 
reduced and for this reason was deemed superior. 

The intrinsic accuracy of the results is +0.2 
km/sec if one considers errors in determining the first 
arrival time of the wavefront at the stations of the 
array and the variations in the azimuth of the wavefront. 

The data from the individual Peace stations were 
combined with station FSJ and averaged after dT/dA com- 
putation and the same operation without averaging was 
applied to data from stations MCC and EDM. 

Figure 15 shows one particular record suite. 
Figure 16 displays the results from VASA in the tele- 
seismic distance range 20° to 100° with the solid curve 
representing the apparent slowness as determined from 
the tables of Jeffreys and Bullen (1958) by cubic spline 
fitting and differentiation. The epicentral distances 


have all been corrected to correspond to earthquakes 
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Figure 15. The vertical component of initial ground 


motion for the Tonga event of July 28, 1970, at 04 h 
47 m and 47.7 s UT as’ recorded by stations oOrecue 
VASA array. The WWVB signal and indicating second 
marks are shown with each seismogram. The magnitude 


was 5.2, and the depth was normal (33 km). 
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Figure 16. The dT/dA observations of initial P waves 
for all events between 20° and 100° as recorded from 
June to August 1970 by VASA. The dT/dA curve is 
obtained from the J-B tables for a normal depth source. 
Squares indicate events from the Aleutians and Japan; 
triangles, events from Tonga and Samoa; plus signs, 
events from South America; circles, events from Asia; 


crosses, special Solomon region events. 
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occurring at normal depth (33 km). As can be seen by 
inspection the majority of events are in agreement with 
the Jeffreys-Bullen (J-B) curve, some of the scatter 
being immediately attributable to the relatively low 
magnitude of some events and the fact that for some 
events part of the array was inoperative. Johnson (1969) 
has previously observed small scale daT/dA anomalies for 
different earthquake regions and such seems to be the 
case for some events shown here. There is, however, a 
striking deviation from the normal in the distance range 
84° to 95° from earthquakes in the Tonga and Samoa 
islands region whereas events from the Solomon islands 
and New Ireland (angular separation from Tonga of 20° at 
the receiver) and Asia and South America show normal 
aT/dA results. Included in the anomalous results are 
one event from the Easter island region and one event 
from Taiwan. The anomalous results have dT/dA values 
that are low (phase velocities therefore will be high 

as shown in figure 17a) with a generally monotonic beha- 
vior with distance. Figures 17b and c show the supporting 
results as determined by the Peace- FSJ array and the 
MCC-EDM station pair. Additionally through the courtesy 
of Dr. D. Davies, an independent analysis was made of 
LASA recordings for 18 events (1966-1969) in the Fiji and 
New Hebrides regions. The phase velocities from this 


study are shown in figure 18 and indicate variations from 
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Figure 17. Observations of phase velocity for Tonga 
earthquakes as recorded by (a) VASA, (b) the FSJ- 
Peace array, and (c) the EDM-MCC station pair. The 
phase velocity curves are for a lower mantle based 
on Birch's (1964) model 2 and its modification (d). 
The solid sections of the phase velocity curves are 
for the geometric rays, and the dashed portions are 


diagramatic for the wave in the high-velocity wave guide. 
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Figure 18. LASA P-wave phase velocities for the 
distance range 91°- 100° and the azimuthal range 
240°- 270° (Fiji and New Hebrides regions). The 
solid phase velocity curve is derived from the J-B 


tables. 


68 


(93a) 


JONVLSIGQ IWHYLNSOIdS 


OOl G6 06 S8 08 


Olc -O¥e 





SJONVY IWHLNWIZV 


S3ILIOOTSA SAVMn-d VSV1 


(93S/WH) ALIDOT3SZA ASVWHd 


Figure 18. 








_3VAW-9 A@AJ 
JAHTUMISA 


pease ewe 


e 





2 
i 
c 






oc: ag 
(.8030) 3OMAT2@IO JANTH3IDNIA 


Bd 


69 


23.9 to 25.7 km/sec at epicentral distances of 91° to 
100°, seven events being low with respect to the LASA 
baseline or J-B derived curve and 11 having slightly 
high phase velocities. Fourteen of the LASA events have 
velocity anomalies that are less than +0.3 km/sec. 
Further, the azimuthal anomalies for VASA are generally 
less than 3° with one being 5.5° and those for LASA are 
generally less than +2° with one as large as 7°. 
Curiously, most of the azimuthal anomalies for VASA 
tend to be clockwise from the geocentric azimuths to 
the sources whereas those for LASA tend to be ina 
counter clockwise sense. More will be said about this 
shortly. 

Upon analysis of the arrival times at individual 
VASA stations one sees a very broad azimuthal effect 
with these stations differing by only -0.2+0.7 sec from 
the J-B tables (corrected for ellipticity and station 
elevation) when all events from epicentral distances of 
20 tO 395° and azimuths 0°) to 360° are considered, The 
analysis of quadrants and individual stations is given 
in table 2. Since many of the events come from Benioff 
zones (see Sorrells et al., 1971) the azimuthal depen- 
dence will be due, in part, to source effects. As can 
be seen individual stations differ by only 0.1 to 0.3 
sec in any single quadrant. For all events into VASA 


the azimuth of arrival was determined independently of 
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Table 2. Differences between observed travel times at VASA 


and the J-B tables for all earthquakes between 30° and 100° 


are shown. 


Renate iE Time Standard 
Station Azimuth Difference | Deviations 
Observations 
(sec) (sec) 


All stations 
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phase velocity and in general showed only slight devia- 
tions from the great circle path. From the examination 
of the small time residuals between individual stations 
it can be stated that structural effects near the 
receiver cannot explain the observed anomalous dT/dA 
values at distances of 85° to 95° from Tonga. Addi- 
tionally there are other events from the Pacific and 
Asia which deviate from the worldwide average values 
for dT/dA as represented by the J-B curve (figure 16). 
Figure 17 contains the observations from the three 
groups of seismograph stations in western Canada with 
the anomalously high phase velocities being observed on 
all three groups in the distance range 84° to 95°. 
Figure 19 shows the epicenters and turning points 
(corrected for depth of source which varies from 30 to 
673 km) for the South Pacific events with the receiver 
systems also shown along with some representative great 
circle paths. As is seen these turning points lie 
within a region 10° in diameter at the end of the 
Hawaiian ridge. Due to the location of the Peace array 
with respect to the TongaSamoa region several normal 
events nee indicated by figure 17b as corresponding 
closely to the Birch 2 curve. These events have their 
epicenters at the southern edge of the Tonga region and 


Since their phase velocities agree with the J-B values 
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Figure 19. Location of earthquakes (plus signs) and 
recording sites (triangles) used in this study. A 

few great circle paths have been drawn to show the 
projection of the seismic rays on the surface of the 
earth. The turning points of rays having anomalously 
high phase velocities are located by crosses. Turning 
points with phase velocities concordant with the J-B 


values are denoted by squares. 
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their turning points (squares in figure 19) appear to 
define the western limit of the anomalous region. In 
addition a single event having a normal phase velocity 
from Fiji into the LASA array has been reported by 
Iyer and Healy (1972) and it has its turning point 
plotted in figure 19. It seems to set a limit on the 
eastern edge of the anomalous region. 

An expansion of this region is shown in figure 20. 
The data used in this diagram are the Tonga-Samoa VASA 
data and the Fiji-New Hebrides LASA data as they include 
both phase velocity and azimuth information. The 
diagram shows both the geocentric turning points of the 
ray (closed circles for LASA and open circles for VASA) 
and the turning points determined by taking the observed 
azimuth of the arriving wavefront and then extrapolating 
back along this ray to the halfway point (indicated by 
the arrowheads). The numbers at the arrowheads refer to 
the actual phase velocity anomaly (in km/sec above the 
J-B values for a particular ray) associated with that 
scattering point. In general when plotting the actual 
turning points those for VASA tend to move to the north- 
west and those for LASA tend to migrate to the southeast, 
thus separating from one another. Further, it must be 
remembered that for these events the distances to LASA 


average 93° whereas for VASA they average 89° and hence 
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Figure 20. Contour diagram of anomalous (observed- 
J-B) phase velocities for earthquakes in the South 
Pacific as recorded by LASA and VASA arrays. The 
tails of the arrows (solid circles for LASA and open 
circles for VASA) represent the geocentric locations 
of the turning points on great circle paths between 
the source and receiver. The arrowheads represent 
the actual turning points based on the observed 
azimuthal deviations. A map of the Hawaiian Islands 
is projected to the level of the turning points which 


are: just above the core-mantle interface. 


76 


HLIWON $334930 


Lv 


‘ONOT J83M $334930 
Srl 


938/ WHO 


ALIDOIZSA 
ISVHd 
AO1 


II WAWH 


INWI 
cs 


®@ ivnvy 


* VOHIN 


SLNIOd ONINUNL GIANISGO — 


SLNIOd SIMWLN39039 VWSVWA eo 
SLNIOd OSIMIN39039 VWSV1 eo 


SLNIOd ONINUNL AVY 





Figure 20. 


ay 


iy S 








EewsEs WOWLH YL 
Ph irre! nak soe y ‘ 


‘og eupthqgnakes in the South 





pon wed WASA areas. (The 
evuawe jos iii ast open 

nA ee 
a ‘icole pathy Setwean 


acrewMecads § 


rous 







MEBL 








: a of tiw testing 
> 
t 


: Aewted wo 4 * : = 

© ire tiovve Es —~ No. tala aaiil : </ ¢ 
bf ° é < 

ee 3 

en, 


DECKEES 





a7 


the rays to LASA will generally bottom as much as 100 km 
below those to VASA. Therefore, thinking three dimen- 
Sionally the two arrays are sampling different regions 
at the ray turning points. The contours are drawn for 
illustrative purposes to help discern the patterns in 
phase velocity distribution: one has a central region 
of high phase velocity (containing points with anomalous 
velocities as high as 3.7 km/sec) centered just north- 
west of the island of Hawaii which becomes more and more 
normal as one moves outward from this region until normal 
or below normal phase velocities are encountered. 

Additional evidence of anomalously high phase 
velocities in this region is furnished by the S-wave 
arrivals. The velocities for the S wave phases were 
determined by the cross-correlation and manual methods 
and also by the coherency method which was discussed in 
Chapter II. The velocities into the VASA array from the 
Tonga islands region are plotted in figure 21 with the 
solid curve representing the J-B derived phase velocities 
for S waves. Behavior of the S phase velocities with 
distance seems to follow the behavior of the P wave 
pares parte for the available common events. 

Three critical areas may be examined for the 
source of the anomalous dT/dA results: (1) source effects, 


in particular the location of the hypocenter on a 
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Figure 21. Observations of S-wave phase velocities. 
The open circles are due to earthquakes in the Tonga 
region (P-wave anomalous) and the crosses are due to 
events from Asia and from the Solomon Islands (P-wave 
normal). The solid phase velocity curve is from the 
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descending plate of lithosphere in the vicinity of the 
Tonga islands (Isacks et al., 1968); (2) effects beneath 
the receiver arrays in western Canada; and (3) the 
vicinity of the turning point of the rays, in this case 
the lower mantle beneath the island of Hawaii. 

Although source or receiver effects cannot be 
completely eliminated as contributors to the anomalous 
phase velocities there are several observations one 
can make which cast doubt on the possibility that these 
effects are very important. The theoretical P wave 
residual field due to a thick slab of relatively high- 
velocity material dipping steeply into the mantle has 
been analysed by Sorrells et al. (1971). These residuals 
vary too slowly with azimuth and the change with distance 
is too low by one or two orders of magnitude to account 
for the anomaly. To generate the required phase velocity 
anomaly one could set up a lateral velocity gradient in 
such a downgoing slab. Since rays into VASA at a dis- 
tance of 90° are separated by only 5 minutes of arc 
they will be merely 0.4 km apart at 100 km from the 
source. In order to produce the anomaly over an array 
150 km on a side one would require that two nearly 
parallel rays travel in two separate velocity structures 
separated by no more than 0.4 km so that the faster 


ray would gain almost 0.5 sec relative to theslower one 
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(A = 90°). Such a rough calculation shows that the 
lateral velocity gradient would have to be greater 
than 1 km/sec per kilometer. This concept is so 
unlikely both from a ray and a wave view point that 
source effects must be dismissed altogether. 

The observed phenomenon can be mimicked by a 
dip of 7° to the southwest (the direction of Tonga) 
on the M discontinuity beneath all stations of the 
arrays but this is at variance with dT/dA from the 
other quadrants as displayed in figures 16 and 17. 
The effect of the dip should be to force increases in 
phase velocities from events in the direction of dip 
and decreases in phase velocities from events in the 
reverse direction. Since many events that appear 
normal velocitywise are not only from the directions 
nearly opposite to the dip but are from regions 
separated from anomalous events by at most 20° in 
azimuth at the receivers, the dipping Moho hypothesis 
must be discarded. Anomalous structures can be placed 
anywhere along the ascending or descending arms of the 
ray but then one requires rather large structures in 
the mesosphere which would lend themselves to previous 
detection. Variations in dT/dA are most easily generated 
at the ray turning point and therefore the anomalous 
velocities will be interpreted on the basis that they 


arise from this region in the lower mantle. 
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re Interpretation 


If, now one assumes that the region about the 
turning point of the ray is the most likely location 
of the dT/dA anomalies observed by both arrays there 
are two basic explanations as to how the anomalies are 
generated. The first, and least plausible, is that 
topography on the core-mantle boundary conspires to 
modify the ray path; the second explanation involves 
lateral and vertical inhomogeneities at the base of the 
mantle. 

tne possibility of undulations on the surface of 
the core was investigated by Vogel (1960) who studied 
hundreds of PcP determinations and found travel time 
anomalies that could be interpreted as arising from 
topographic variations of about +100 km on the surface 
of the core. However, such variations are highly suscep- 
tible to uncertainties in source depth and location and 
also to inhomogeneities in the crust and upper mantle. 
There have been several attempts at explaining that the 
observed gravitational field of the earth is due in part 
to such undulations of the core-mantle interface. Hide 
and Horai (1968) have produced hypothetical contour maps 
based on a spherical harmonic analysis of the gravity 
Paella fortdecrec "i = "4y ns 6, and*n S 8 being generated 


by bumps of amplitudes +4, +15, and +66 km respectively, 
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assuming a density increase of 5.57 to 9.74 g/cm? 


across the boundary. 

However, it is now known that the radius of the 
core has been determined to within about 10 km. 
Jeffreys (1939), Taggart and Engdahl (1968), Johnson 
(1969), and Buchbinder (1971) obtained core radii of 
3473+4 km, 3477 km, 348142 km, and 3479+2 km respec- 
tively by using PcP phase observations. Jeffreys' 
value has been used in this study. On the basis of ScS 
observations Hales and Roberts (1970) have determined 
the radius to be 3486+5 km. A very sensitive method 
of measuring not only the radius of the core but also 
determining the existence of any large scale variations 
in the shape of the boundary is afforded by multiple 
reflections such as P3kP, P5kP, P7kP, and so on, where 
the numeral refers to the number of internal reflections 
suffered by the compressional wave inside the core. 
From 20 observations of the P7kP phase Buchbinder (1972) 
has found a model that agrees to within 1 sec of the 
data. Therefore, it is apparent that observations of 
Gata associated with the core tend to rule out the 
ahi sebace of bumps of more than a few kilometers in 
amplitude in this boundary. In order to produce the 
observed phase velocity anomaly ray tracing shows that 
one would require a depression in the core of nearly 


125 km with an increased mantle velocity gradient. 
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However, this would clearly lead to a negative gravity 
anomaly which is just the opposite of the observed 
field and thus eliminates this model. 

Davies and Sheppard (1972) have done an extensive 
study of dT/dA at many distances and azimuths. Their 
study, using the large aperture seismic array in Montana, 
indicates a rapid change in dT/dA in the vicinity of the 
core-mantle boundary near Hawaii, the Galapagos islands, 
and Iceland. Julian and Sengupta (1973) using travel 
time evidence have also found anomalous structure 
beneath Hawaii in this depth range. Some of their 
rays seem to sample high velocities and some sample low 
velocities although an explicit spatial distribution with 
respect to the islands has not been given. Their sampling 
also indicates that this anomalous structure extends 
vertically several hundred kilometers. 

The results of this work lead to the interpreta- 
tion that there exists a lateral and vertical hetero- 
geneity near the core-mantle boundary the surface projec- 
tion of which lies northeast of the island of Hawaii and 
is approximately 5° in diameter, The high velocity core 
of this heterogeneity seems to be surrounded by a normal 
to low velocity aureole which may constitute the plume 
material as it is heated by the central, possibly rigid, 


structure causing it to flow upward. The interpretation 
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assigned to velocities may be justified as a first 
approximation since the actual azimuthal deviations 
suffered by the ray from the great circle path are 

small (<5°). The velocity model used as a standard 

for this study was Wiggins' (1968) self-consistent 
velocity and density distribution which was based on 
model 2 of Birch (1964). A continuous velocity function 
within the earth's mantle was generated by fitting the 
discrete values with a cubic spline interpolation for- 
mula. The smoothing spline routine (ICSSMU) of the 
International Mathematical and Statistical Libraries 
based on the work of Reinsch (1967, 1971) was used in 
order to avoid numerical instabilities. The geometrical 
ray integrals (Bullen, 1963) were calculated using a 
computer program based on the work of Chapman (1971). 
Upon examination of the three Tonga-Samoa data sets 
(figures l/7a, b, cc) it’ is clear that the experimental 
phase velocities start to deviate from Birch's model 2 
near 84° and continue in this manner smoothly to 91°, 
beyond which there is a marked increase in the observed 
velocities and also in their scatter. Birch's model 2 
has been er iere to include a higher velocity gradient 
between depths of 2371 km and 2861 km (figure 17d) so 
that a fit to the data in the range 84° to 91° was made. 


The higher phase velocities for distances beyond 91° 






jeri? 6 8s betiiteut sd yam esitrcolevy of boaptees 

evoidetveb Leddomine tso9oe ode Sonte noltemtxorggs 4 

e15 ieq sioris Jserp sft Mega yer on? yt bevetiwe 
brsbeste 8 aa beer Lebom ytiselev sdt  .(*e@r) Clam 
tastelecoo-tise (808!) “eabpplW es¥ ybode eid? 102 
co besed sew dotdw noltodiisgeth yeiensb bas ysisoley 
nokzonut ysatsoolsey evountines A’ O¢RRRE) desk to f fesoe 
oct piisi2i yd Bes ~~ asw sitnem a'trse eft aldgiw 
-1o1 nuoiteloyxednt satiea tidus s diiw eenlay steto8rs 
ofd to (uMeeDL) entsaon ondlee padvtoone eft stom 
2otterdid Leoidyidss® bas Isobieeenset Canoiienresnt 


heave asw ({0eL ,T8CL) foentiedt Yo tyow edd ado Beged 


(sovincardeap oT .eeirvitcdsdenn Iscigeha Sinve Oo tebh20 
‘tev beveluoles stew (£081 ,aellug) eiempegal vax - : 
(TOL) mamged9 Yo Wav sie do beaded mexpoay seeeqneo ; 

atee at6b somsd-apne? eerdt sa te dotgeaigee aeg7 
lesnemiteaqxs of% Feds vests et 22 (5 Va sees asxupit) 
¢ febom a'fortle moxt ofetveh of S¥saa aslzksoisv easag 
,°1@ of yttidoome xsftnam eina ok Sdnigae> Bae “$8 7e0n - 
bevisedo oft ai azes roni betism 6 ab oxsdd fobiw baoyed © _ 


$ [bom e‘dowis 1933892 1f0o09 nib o8lf6 She eetstoolev ’ 
Serpe, apatey reigns Camndieaeaalll 


ee: Ss nah eat t: 


vais, 





86 


required a sharp increase in mantle velocity to at least 
14.7 km/sec between 2861 km and 2898 km. A smooth 
increase in velocity to 14.7 km/sec at the core-mantle 
boundary is unacceptable since the shadow would occur 
at too small an angle and the decay rate of the signal 
would be increased (Phinney and Alexander, 1969; 
Chapman and Phinney, 1970). The observations beyond 
about 89° would not be of a direct P wave but of a 
trapped interface, or 'head' wave type. Evidence of 
this is visible in figure 22 which compares an anomalous 
event with a normal event both from the South Pacific 
area (Tonga and Solomon islands). The anomalous event 
exhibits a distinct low frequency precursor which is 
very much like the head wave precursors generated by 
Waddington (1973) who has used the Cagniard-de Hoop 
method to produce synthetic seismograms for this distance 
range. Of all his velocity-depth models one similar to 
that shown in figure 17d gives rise to seismograms that 
are very much alike in appearance to the observed wave- 
forms. The very high phase velocities shown in figures 
1l7b and c for distances less than 88° may be due to 
pac cetian (dashed lines) from the high velocity boundary 
layer. 

Both the azimuth-velocity anomaly distribution 
of figure 20 and the scatter of the phase velocity 


diagrams of figure 17 indicate that the actual structure 
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Figure 22. The wave forms of two events, one from the 
Solomon Islands having a normal phase velocity and one 
from the Tonga Islands having a high phase velocity. 


The events are each shown for three stations of VASA. 


SOLOMON ISLANDS EVENT: NORMAL 
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of the anomalous region is quite complicated and may 
contain further layering in a vertical sense along 
with lateral velocity variations as evidenced by the 
transition from high velocities to low velocities as 
the structure is traversed. The anomaly is almost 
certainly no more than 300 km wide and extends upward 
about 300 to 400 km (from the phase velocity diagrams 
and the maximum depth of penetration of these rays), 
although it is very difficult to find an earthquake- 
receiver system to adequately determine the upper limit 
on the vertical dimension. 

Figure 23 shows a cross section of the earth 
with the anomaly sketched in at its proposed location 
at the base of the mantle. The resultant plume is 
also shown rising through the mantle to the surface 
where the volcanic islands are generated on the moving 
Pacific plate. In this diagram are shown a number of 
seismic rays emanating from what would be a source in 


the Tonga region. 
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Figure 23. Sketch of the anomaly at the base of the 
mantle beneath the Hawaiian linear volcanic chain 
showing the location with respect to the rays passing 


through this region. 
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CHAPTER IV 


THE HAWAIIAN LINEAR VOLCANIC CHAIN AND ITS POSSIBLE ORIGIN 


7 PLncrOoaucLiLon 


Upon examining the earth's island systems it 
becomes apparent that the Hawaiian Islands taken with 
the Emperor Seamounts are representative of several 
linear volcanic chains throughout the world. For 
example one notes the presence of the Austral-Marshall 
and Tuamoto-Line chains in the Pacific, the Rio Grande 
and Walvis ridges emanating from Tristan da Cunha and 
Gough Island in the South Atlantic, and various volcanic 
provinces in Africa and Europe. In the light of recent 
mantle hot spot and plume hypotheses it is rather 
desirable to examine the lower mantle beneath these 
regions in order to determine the nature and location 
of any features there. This, combined with crust and 
upper mantle data,will provide a picture of the nature 
of the earth over a broad range of depths which is 
essential for any investigation of the accuracy of the 
plume hypotheses, their possible worldwide occurrence, 
and their relationship to the earth's plate tectonic 


system. 
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2. Surface Description of the Hawaiian Linear Volcanic 

Dana (1890) who conducted a rather extensive and 
detailed survey of the Hawaiian Islands during the middle 
of the 19th Century recognized the essential volcanic 
nature of the Hawaiian system and described in some detail 
the processes of volcanism in this area. He recognized 
also that the chain is bounded on both sides, at least in 
the vicinity of Oahu and Molokai, by a more or less con- 
tinuous trough and that this was probably due to post 
volcanic subsidence. 

Jackson et al. (1972) give a very lucid and 
extended description of the Hawaiian Archipelago. They 
view the linear chain as being a more or less continuous 
sequence of shield volcanoes stretching across the 
Cretaceous sea floor from the island of Hawaii to Midway 
Atoll and then bending toward the northeast in a line 
forming the Emperor Seamount chain terminating in the 
Aleutian trench. Jackson et al. describe the 3500 km 
ridge as being composed of shield volcanoes rising as 
much as 9 km above the sea floor and 5 km above sea level 
having diameters as great as 120 km, all capping the 
Hawaiian Ridge, a well defined topographic high on the 
ocean floor bordered on both sides by a moat up to 700 m 
deep. The individual shields are separated from one 


another by an average distance of 75 km and began with 
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the copious eruption of tholeiite (the gentle, free- 
flowing nature of eruption has led volcanologists to 
characterize all such eruptions as being of the 
"Hawaiian type") and was followed by a capping of a 
much smaller amount of alkalic lava. They point out 
that projected linear magnetic anomalies cross the 
Hawaiian chain at angles varying from 30° to 60° and 
cross the Emperor chain at 80° to 100° showing no change 
in trend from one side of the chains to the other. The 
available evidence, they note, indicates that the orien- 
tation of the chains has not been influenced by the 


structure of Cretaceous floor. 


3. Geochronology of the Hawaiian-Emperor Chain 


The first systematic investigation of the ages 
of the Hawaiian Islands was done by Dana (1890) who found 
evidence of a geomorphological nature that the chain 
became younger from northwest to southeast by observing 
the progressive state of ruination of the volcanoes in 
the reverse direction. 

This trend has been quantitatively substantiated 
by McDougall (1964) who has examined samples by K-Ar 
dating from Kauai to Hawaii and found in general that 
mean ages for lava deposits decrease toward the south- 


east (Table 3). 
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3. Ages of the Hawaiian Islands as a function of 


distance from the active volcano Kiluea on Hawaii. 
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The trend of increasing age as one moves from 
Hawaii to Kauai and Midway is supported by Funkhouser 
et al. (1968) who has also employed K-Ar dating to 
determine an age of 7.5 my for Nihoa Island and 11.3 my 
for Necker Island. Paleontological investigations by 
examination of drill cores on Midway Atoll has established 
that the island was formed in pre-lower Miocene times 
having an age greater than 25 my (Ladd et al., 1967). 
Recent microfossil evidence taken from the most northerly 
feature of the Emperor Seamount chain, Meiji Seamount, 
indicates’ an’ age of°about 70 my (Schollret al.,'1971; 
Berggren, 1972). Clague and Dalrymple (1973) have 
examined material dredged from the top of Koko Seamount 
lying 300 km north of the Hawaiian-Emperor bend, the 
mean age being determined to be about 46.4 my. They 
further describe this seamount as having the typical 
smooth form of the Hawaiian volcanoes probably being 
formed by the coalescing of several shields. 

Thus while extremely sparse as yet, the existing 
age data from the seamount chain tend to support the 
trend established by the much more extensively sampled 


Hawaiian chain. 


4. Subsurface Data of Hawaii 


Gravity data has been examined by Gaposchkin and 


Lambeck (1971) who have used satellite data to produce 
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a graph of the global free-air gravity anomalies with 
respect).to. the, best fitting ellipsoid.. The gravity 
field to the sixteenth degree shows the Hawaiian region 
within a broad positive anomaly. Kaula (1972) has also 
noted that the broad positive gravity anomaly around 
Hawaii is the only one of its type not associated with 
a spreading center. He further indicates that there are 
areas having this type of broad positive free-air or 
isostatic anomaly which seem to be connected with rise 
formations: Yellowstone, Galapagos, Easter, Eltanin, 
Balleny, Prince Edward, Bouvet, Iceland, Azores, Afar, 
Kamchatka, among others. 

Pekeris (1935) considers two opposing contribu- 
tions to gravity anomalies due to the rising part of a 
convection cell. There is a positive contribution to 
gravity due mainly to the rise of the upper surface of 
the crust and there is a negative contribution to gravity 
on account of the lower density of the matter being 
pushed up caused by the elevated temperature of the 
rising stream. He points out, however, that the positive 
contribution is the larger. McKenzie (1967) states that 
the long-wavelength harmonics of the gravity field as 
determined by satellites are not related to shallow 
crustal.or, lithosphere, effects, but.are probably related 


to flow patterns within the deeper mantle, and Runcorn 
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(1965) argues that departures from hydrostatic equilibrium 
as evidenced by satellite observations of the geopotential 
are due to flow within the mantle. 

Doell and Cox (1972) have used the idea that the 
nondipole magnetic field generated within the earth's 
core may be used to investigate lateral and vertical 
inhomogeneities at the core-mantle boundary. These 
variations in the lower mantle will generate differences 
in the secular changes of the field and the strength of 
the nondipole field at points of observation on the 
eantipes sunface..\,,Doell iand Cox ,point,.out,.that it is not 
attenuation in the lower mantle that gives rise to these 
variations in the longer period geomagnetic spectrum but 
rather, lateral changes in the source function. The 
nature of the source function can be controlled then 
by properties of the lower mantle by interface coupling 
with the core through the core-mantle boundary. From 
paleomagnetic and geomagnetic observatory data they have 
found a subdued secular variation and a subdued nondipole 
field in the central Pacific consistent with a pronounced 
attenuation of the geomagnetic spectrum in the period 
range 200 to 2000 years. This implies coupling exists 
between a lateral heterogeneity in the lower mantle and 
the core resulting in partial extinction of the non- 
dipole fluctuations in this period range. Limited know- 


ledge of the manner of magnetic field generation in the 
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core prevents inversion of the data to find a suitable 
model and also prevents more detailed analysis of other 
more complex regions. 

In a series of papers (Alexander and Phinney, 
1966; Phinney and Alexander, 1966, 1969) Phinney and 
Alexander pursued a study of the core-mantle boundary 
region by the examination of the spectra of long period 
diffracted P waves (since these are the first to arrive 
in the shadow region). The source radiation effects were 
removed by taking the ratio of the spectra obtained from 
the vertical ground motion at two stations at least 20° 
apart and on the same azimuth to the source. The results 
for waves diffracted at the core-mantle boundary beneath 
Hawaii showed a broad and pronounced absorption peak in 
the frequency range 0.03 Hz to 0.04 Hz with smaller 
peaks being observed at still lower frequencies. Phinney 
and Alexander considered the kinematics of the deep 
shadow diffracted wave which travels on the boundary as 
an interface wave offering sensitivity to lower mantle 
structure. They attempted to interpret radial and lateral 
variations observed in terms of layered lower mantle 
models 30 to 160 km thick having weak absorption peaks 
in the decay spectrum and noted that for paths where the 
decay is pronounced this corresponds to a positive 
velocity gradient of 0.2 km/sec per 100 km. At short 


periods they suggested that this gradient may be even 
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greater. For rays diffracted under the Atlantic Ocean 
off Africa they found no such decay pattern. Although 
Phinney and Alexander did not produce a particular 

model to fit their Hawaiian observations they did 
pioneer the study of lateral and vertical heterogenei- 
ties near the core-mantle boundary. Mitchell and 
Helmberger (1973) using the amplitude ratio of long 
period SH and core reflected SH, ScSH/SH (taking the 
ratio to minimize source and sub-receiver effects), have 
produced a layered model at the base of the mantle. 
Comparison of the observed with the synthetic seis- 
mograms computed by the Cagniard=de Hoop method indicates 
that neither the Jeffreys-Bullen nor a negative velocity 
gradient at the base of the mantle satisfies the obser- 
vations. Rather, a positive velocity increase of between 
0.3 and 0.7 km/sec over a depth range of 40 km to 100 km 
is deemed necessary. Further restrictions on this model 
are determined by considering the differential times, 
ScS-S, for both radial and transverse components, the 
latter being smaller by an amount explainable by a high 
velocity region at the base of the mantle 20 to 70 km 
in thickness having a velocity increase of between 0.3 
and 0.5 km/sec. It should be emphasized that the long 
periods of the shear waves tended to mask lateral 
heterogeneities (which was desired by the authors). 

This study, nevertheless, does indicate the existence of 


positive velocity gradients in this depth range. 
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5. Origins of the Hawaiian Linear Volcanic Chain 


Several explanations of the genesis of the Hawaiian 
linear volcanic chain have been put forward. Wilson 
(1963) proposed that a source of lava rising from a 
relatively stagnant center of a plate-driving convection 
cell in the mantle would puncture the lithosphere depo- 
Siting volcano producing material in the crust as the 
plate was pushed along thus forming a line of successively 
younger volcanoes in a direction opposite to the sense of 
plate motion. 

Dana (1890) using the old supposition that vol- 
canoes arise over the intersection of two faults states 
that the Hawaiian Islands were formed by two parallel 
northwest-trending "fundamental rifts" along the course 
of the island chain with individual volcanoes being 
generated where these rifts were crossed by northeast- 
trending fractures. This view was supported by MacDonald 
and Abbott (1970) whose explanation differs little from 
Dana's, the orientation of the supposed rift being 
moderately disrupted by slippage along the intersecting 
Molokai fault. McDougall (1971) and Green (1971) modify 
this explanation to a progressive intrusion of magma 
along lines of structural weakness. 

Morgan (1972) gives an explanation similar to 


that of Wilson (1963) except that his source or "hot spot" 
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is placed in the lower mantle where it generates a 

plume of material rising to the earth's surface. 

Jackson et al. (1972) have reviewed evidence that 

tremor swarms associated with magmatic movements and 
eruptions are highly localized being confined to a 
cylindrical region 25 km in diameter and 25 to 60 km 
beneath the volcano which is compatible with the shield 
characteristics of these volcanoes. Hence these 
extremely localized sources tend to rule out the rather 
broad upwelling of material inherent in the rift hypo- 
theses. Ryall and Bennett (1968) have examined the 
crustal structure of Hawaii in great detail finding it 
complex with the presence of transcrustal fracture zones 
related to uplift and volcanism. They find, however, 
"no evidence for large-scale transcrustal faulting as 

a controlling factor for volcanism along the Hawaiian 
Ridge' (Rvall and bennett, 19668, .p.4561) 2 -Hiil (1969) 
has done a large amount of seismic refraction work 
involving the crust and upper mantle underlying the 
island of Hawaii and has noted that the upper mantle 

P wave velocity probably decreases beneath Kilauea. 
Schilling (1973) has indicated that these plume regions 
such as Iceland and Afar show similar low velocities and 
that this is to be expected for density-deficient material 


upwelling from deeper in the mantle. 
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6. Worldwide Nature of Hot Spots and Controversies 


Since Wilson (1963) proposed that the genesis of 
the various volcanic chains is due to the motion of 
lithospheric plates over hot spots fixed within the 
mantle with respect to each other, Morgan (1972) has 
provided a list of 20 probable hot spot locations 
(figure 24) each being associated with a mantle convec- 
tion plume. He has noted the worldwide correlation of 
Hawaii type positive gravity anomalies with the postu- 
lated hot spot locations and states that they are 
symptomatic of rising currents in the mantle: the less 
dense material of the ascending plume produces a broad 
negative anomaly which is offset since the surface of 
the crust is deformed upward causing excess mass to be 
closer to a satellite passing overhead making the net 
anomaly positive. 

A test of the hypothesized fixed nature of the 
plumes is whether the directions of the volcanic traces 
coincide with the direction of instantaneous motion of 
the plate and whether the age determinations support 
the accepted plate speeds over the past several million 
years. Such an investigation has been carried out by 
Minster, Jordan, Molnar, and Haines (to be published) 
who have numerically modelled the instantaneous plate 


tectonics. Their observations of the Pacific, North 
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Figure 24. Worldwide plume locations (Morgan, 1972). 
The abbreviations are as follows: 

Leek - Afar 

2. Am - Amsterdam 

3. Asc - Ascension 

Az.) Az w= JAzZoOres 

5s 1Bal = )Bableny 


6. Bou - Bouvet 


LE ae - Canary Islands 
Some aks - Comores Islands 
Oe - Hawaii 


lO “wee. = lee land 

ele oss - Mount Kenya 
12. McD - McDonald 

13. PE - Prince Edward 
14. Re - Reunion 

T5u" sStewhy— Steshetena 


16. Tr - Tristan de Cunha - Gough 


The names in bold face print identify the various tectonic 


plates. 
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American, South American, African, Eurasian, Antarctic, 
Nazca, and Cocos plates indicate that the fixed plume 
assumption is indeed in accord with the vast perpon- 
derance of data available. In the case of the swiftly 
moving Pacific plate, for example, their model predicts 
a present rate of motion over the Hawaiian plume of 8.9 
cm/yr at an azimuth of N 66°W compared with their quoted 
observed azimuth of N 64°W and rates of 10.9 cm/yr 
(BoLcdi, 23/2), 10 Gmyvr (MCboucall, 1964), 9 cih/yr 
(Malahoff and Woollard, 1970) as determined by radio- 
metric dating. 

Duncan et al. (1972) have found that assuming the 
Eurasian plumes are fixed one finds a large discrepancy 
between the plume polar wander curve and the actual paleo 
pole positions. They state that this may imply either 
shifting of the entire lithosphere or a movement of the 
plumes. McElhinny (1973) has shown that the sum of the 
vectors (on an equatorial plane) of plate motions rela- 
tive to the earth's rotation indicates that the litho- 
sphere has not moved as a whole in Tertiary times. The 
results of Grommé and Vine (1972) who have studied paleo 
latitude measurements show the same latitude for Midway 
Atoll as for present day Hawaii within experimental 
error. Thus, if no true polar wandering has occurred 


then this plume, and probably the others, must be fixed. 
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The paradox is resolved, however, since Burke 
et al. (1973) have shown that Duncan's interpretation 
of the Eurasian traces upon which his results are based 
is in error: they describe two plume traces near 
Iceland-Myvath to Faeroes and Heckla to Rockall Bank, 
whereas Duncan has incorporated these into a single 
somewhat meandering trace leading to Scotland; also 
Burke et al. regard the Central European volcanic 
province as being composed of two traces both north- 
south trending in disagreement with Duncan's earlier 
assumption of one trace running east-west. 

Morgan (1972) has given calculations of the ver- 
tical and transverse stresses exerted on a plate due 
to a rising plume. The gravity anomaly created by a 
plume with a diameter of 150 km is of the order of 20 
mgal falling to zero at a distance of 1000 km along 
the surface which is in accord with actual observations 
(Kaula, 1973). The shear stress produced on the under- 
side of the plate by such a plume (rising at a rate of 
2 m/yr) would be sufficient to drive the plate along at 
about 5 cm/yr. By suitable choice of viscosity and 
plume flOw rate the total force exerted on the plate can 
be greater than the resistive drag on the plate bottom 
and along the edges. The key estimate is the volume of 
flow of the plume. Conceivably, plumes situated 
at rise formations may then represent the driving mech- 


anism of the earth's plate tectonic system. 
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7. Conclusion 


From the evidence presented in Chapters III and 
IV it is concluded that the lower mantle does exhibit 
lateral and vertical inhomogeneities. An anomalous 
region possessing a high velocity central zone with 
a normal to low velocity aureole appears to be present 
just above the core-mantle boundary underneath the island 
Of §Hawaaie (Ligure, 20; Chapter, LIL) ».-The»surface, projec- 
tion of this region, in fact, is located in an area where 
Jackson et al. (1972) have placed the position of the 
present Hawaiian plume from geological considerations. 
The time of formation of this heterogeneity must have 
been early in the history of the earth's development 
when the core was being differentiated. 

To the present time there has been little satis- 
factory evidence supporting any of the competing hypo- 
theses for the origin of the Hawaiian linear volcanic 
chain as outlined in this chapter. Neither a fundamen- 
Sales NOL a priamary source tor the Javas has, been 
established. The asthenosphere is an unlikely source 
in light of plate tectonics and Wilson's (1963) hypothesis 
and from seismicity studies on the depth of earthquakes 
and seismic observations of mantle structure there is 
no indication of a source in the mesosphere. Hence it 


is indeed possible that a plume, whose surfacing point 
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lies just northeast of the island of Hawaii, does exist 
and that its source lies in the anomalous high velocity 
region at the base of the mantle. This picture, although 
not unequivocally proven as far as a direct connection 
between the region at the base of the mantle and the 
surface is concerned, is not at variance with the dT/dA 
observations. It should be noted that Shaw and Jackson 
(1973) prefer an alternate theory for Hawaiian type 
melting spots. They speculate that shear melting caused 
by plate motion leads to the formation of a dense residuum 
which descends to the core-mantle boundary forming a dense 


high velocity region. 
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THIS PROGRAM IS A GENERAL VELOCITY —- AZIMUTH SPECTRAL 
ANALYSIS PROGRAM PRODUCED BY P. Re. GUTOWSKI AT THE U. 

OF ALBERTA. 

ese LATEST UPDATE) ——) MAR L197 73) os 

THE ORIGINAL VELOCITY SPECTRAL ANALYZER WAS GIVEN BY 

De DAVIES ET AlLes NATURE 1970+ AND FEATURED BEAM FORMING 
FOR EACH VELOCITY (DELAY AND SUMMATION) AND THEN FOUND 

THE POWER IN 1 SEC OF BEAM IN 1 SEC INCREMENTS DOWN THE 
RECORD. THE RESULTING VELOCITY VS. TIME MATRIX WAS THEN 
CONTOUREO ANC PLOTTED. THIS METHOD HAS SINCE BEEN SHCWN 
TO BE RATHER DEPENDENT ON AMPLITUDE VARIATIONS ACROSS THE 
ARRAY WITH DISTINCT SIDE LOBES APPEARING EVEN FOR LARGE 
NUMBERS OF SENSORS. THE VELOCITY — AZIMUTH SPECTRAL NETHOD 
(COVESPA) COES NOT SUFFER FROM THIS CRAWBACK AS IT INVOLVES 
A COHERENCE FUNCTION GENERATED BY CROSS MULTIPLICATION OF 
TRACES IN COMBINATIGNS OF TWC NOT THE SIMPLE SUMMATION OF 
THE SENSORS I.E. FOR S SENSCRS BEAMFORM WOULD LEAD TO THE 
SUM OF S TRACES WHEREAS COVESPA WOULD CROSS MULTIPLY 10 
TIMES (5 C 2) AND THUS MAKES MUCH MORE EFFICIENT USE OF 
THE AVATLABLE DATA REDUNDANCY. IN ADDITIGN COVESPA SWEEPS 
NOT ONLY THROUGH VELOCITY AND TIME. BUT ALSO THROUGH 
AZIMUTH THUS EMPLOYING THE MAXIMUM INFORMATION INHERENT 

IN ARRAY DATA. FOR A MODERATE NUMBER OF SENSORS THEREFOREs 


THERE WILL BE VERY LITTLE SIOE LCBE PROBLEM. 
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INPUT 
CARD1 THRU 14 -- ARAY GEOMETRY (OISTANCES AND AZIMUTHS) 
CARD15 -- FILE NUMBER OF DATA ON TAPE 
CHANNEL O —- VERTICAL» 1 - NSs 2 — EW 
NUMBER OF STATIONS FOR THIS EVENT 
LENGTH OF TIMECSEC) INDEX OF CORRELATIGN MATRIX 
STATION NAMES IN LITERAL FORMAT 
FORMAT. eo 415 » SA4 


NOTE: NFILE=999 TERMINATES EXECUTION OF PROGRAM 


CARD 16 START TIME OF CORRELATION GIVEN IN MIN. 
PROGRAM SEARCHES + OR —- 10 DEG ON EITHER SIDE 
OF THE STARTING AZIMUTH (GECOCENTRIC AZIMUTH) 


THE START TIME OF ALL THE RECORDS ARE IN SECONDS 


CARD17 LOWER SLOWNESS VALUE IN SEC/DEG 
HIGHER SLOWNESS VALUE IN SEC/DEG 
SLOWNESS INCREMENT E.Ge 0220 SEC/DEG 
NOTE: THESE SHOULD BE CHOSEN SO THAT THERE ARE 
LESS THAN 32 SLCWNESSES 
CORRELATION LIMIT BELOW WHICH MAXIMA IGNORED 


FORMAT 4F10.0 
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THIS PROGRAM WILL NOW SEARCH OVER THE SPECIFIED VELOCITY» 
TIME AND AZIMUTH RANGES GENERATING A SERIES OF COVESPAGRAMS 
FOR EACH AZIMUTH PRINTING THESE CUT AND PICKING SUCCESSIVE 
MAXIMA OVER A 5S SEC TIME WINDOW INCREMENTED AT 2 SEC 
INTERVALS THUS PRODUCING A SERIES OF VELOCITY» TIMEs AND 


AZIMUTH ESTIMATES BY PARABCLA FITTING 


WEAK TE MA I We I A A SK a a a a a aK a RK OK i a a a aK a ee aK KO a OK aK a OK cK a 2 a I a OK I IK IK 


DIMENSION AZ(7s7)2D(7s7) 9V( 592000) sCC(10510) 

DIMENSION CCOR(120231421 )sDELTA(10s10) sAZM(7) »TIMES(S) 
DIMENSION XVEE(31) 

DIMENSION AT(21) 

DIMENSION IA(7) oL010) »M(10)-sITO(S) 

INTEGER*4 TD(10910) 

REAL*4 STATNS(5) 

REAL*4 KEYTBL(7)/* DEL*%s* EDM's® MAR*%s*® PIN's® RM1%y 
1* RM2%,* TROPS 

INTEGER*2 IDAT(8192) 


NFLST=1 


READ ARRAY DOISTACES FROM STATION TO STATION 


READ(S5,101) O 


READ ARRAY AZIMUTHS FORWARD AND BACK 
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100 


102 


105 


402 
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READ(Ss101) AZ 


FORMAT(7F10.6) 


READ NFILE sCHANNEL sNUMBEROF SENSORS»CORRELATION LENGTH» 


ANDO STATION NAMES 


READ(S5+s100) NFILEsLLsNSTATec I TSPANsSTATNS 
FORMAT(4I5 0S5A4) 


IF(NFILE .EQ. 999) GO TO 9S9 


READ START TIMEs AZIMUTHs AND TRACE START TIMES 


READ(5e102) TSTRTsAZMTHs TIMES 


FORMAT(7F10.0) 


READ LOWER SLOWNESS LIMITs UPPER LIMIT+ CORREL. THRESHOLD 


READ(5s105) VL eVUsVDIV»sCORLIM 
FORMAT(4F10.0) 


ITF GNFILE NE. NFLST) GO TO 401 


IDENTIFY STATIONS 


DOATT TISLINSTAT 

DO 72 J=1.7 

IF (STATNS(I) .NE. KEYTBL(J)) GO TO 72 
IACI) =J 


GO TO 71 
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72 CONTINUE 


71 CONTINUE 


GENERATE ALL POSSIBLE COMBINATIONS OF TWO STATIONS 


AN=NSTAT 


NUM=AN* (AN-1.)72. 


DO 19 I=1»sNUM 
K=K+1 
TRUCKS eCE. NSTAT) GO TO. 21 
J=Jt+1 
K=J+1 
11 LCI )=J 


19 M(I)=K 


FIND AMOUNT OF DATA TO TAKE FROM TAPE ANO WHERE TO START 


TMAX=TIMES (1 ) 
IMAX=1 
DO 20 f=2eNSTAT 
IF (TMAX «GT. TIMES(I)) GO TO 20 
TMAX=TIMES(T} 
IMAX=I 
20 CONTINUE 
DO 200 I=1eNSTAT 


TDEL=TMAX-—TIMES(I) 


cul 
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200 


814 


818 
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IDELAY=12,.5*TDEL 

ITO (I )=4*I DELAY 
BLK=TSTRT/2, 73067 
MBLK=BLK 

XBLK=MBLK 

TREM=2. 73067%*( BLK-XBLK) 
IREM=4*x IF IX ( TREM*750.) 
ITSPAN=I TSPAN#+40 
TSPAN=ITSPAN 
ISPAN=12.5*TSPAN 


WRITE(65814) ITO 


FORMAT(5Xe "ITD. oe § 9 5110) 
WRITE (62818) BLK eXBLKsTREMsIREMs ISPAN 
FORMAT (SX eee *s SF10 0423110) 


DO 210 [=1eNSTAT 


READ STATION HEADER FOR THIS EVENT 


REAC(&) NFe TEVNT sSYSeSTAT 


RE AC(8e END=220) 


SKIP TO STATION DATA FOR THIS EVENT 


CALL SKIP( 0 oMBLKe 8) 
ICOR=0 
READ(8) IDAT 


WRITE(6820) I 
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820 FORMAT(SXs *STATNee el * 0110) 


PICK CFF DATA FOR EACH STATICN SCG THAT ALL TRACES ARE 


ALIGNED IN TIME 


DO 230 J=1+ISPAN 

14=4%J 
K=14+ITREM+I TOC 1)-3+LL+ICOR 
V(IeJ)=IOAT(K) 

IF(K .LT. 81894+LL) GO TO 230 
WRITE (60820) K 

WRITE (65820) 14 

READ(8) IDAT 
ICOR=-IREM-ITO(I)-14 


230 CONTINUE 


REMOVE OC LEVEL 


DC=0.0 
DO 233 J=1+*ISPAN 
233 DC=DCtHV(IsJ) 
DC=DC/FLOATC(ISPAN) 
DO 234 J=1+sISPAN 
234 VC IeJ)=VCIsJ)-DC 
WRITE(69900) NFeIEVNT» STAT 
900 FORMAT(S5SXs*PARTIAL TAPE HEACER *9215sA10) 
WRITE(65821) DC 


821 FORMAT(S5Xe *0C. ee" 9 F102 3) 
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CALL SKIP(1208) 
NFULST=NFLST+2 
ITSPAN=ITSPAN-20 
RAD=3.14159/180. 
IVL=VL*¥100. + 0.2 
IVU=VU*100. + 0.2 
IVDIV=VDIV*¥100. + 0.2 
AZML=AZMTH-10. 

AZMU= AZMTH+4+10. 
IAZML=AZML*10. 
TAZMU=AZMU*10. 

1Z=0 

DO 107 ITAZM=IAZML eI AZMU510 
IZ=I1Z+1 

ATC IZ )=FLOATCIAZM)/10. 
A=AT(IZ) 

AZM(2)=A 

AZM (1)=At+0.17850 

AZM (3)=A-0 .648 
AZM(4)=A-0.08083 

AZM (5) =A-1.63950 

AZM (6 )=A-1.54783 


AZM(7)=A40. 22184 


CALCULATE DISTANCES FROM STATION TO STATION THE WAVEFRONT 


HAS TO TRAVEL FOR EACH AZIMUTH 


DO 250 I=1+NUM 
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DELTACLC IT) eMC IT) V=OCTACL( I) + TACMCT) ) )*COS(CRAD* (AZMCIA(LCI 
TL) I-AZCITACL (CT) ) oTACMOT)))) 


250 CONTINUE 


FOR THIS AZIMUTH CALCULATE COVESPAGRAM 


Iv=0 

DO 306 IVEL=IVL>sIVUs,IVDIV 
IT=0 

IV=IV+l1 
VEX=FLOATC(CIVEL)/100. 
IVLIM=IV 

VEE=111.2/VEX 


XVEECIV)=VEX 


CALCULATE DELAY TIMES FOR EACH 2 STATION COMBINATION 


FOR EACH VELOCITY 


DG 301 I=1.sNUM 


301 TOC LCI) oMCI) )=12.5*(DELTA(L( I) »M(I))/VEE) 


STEP DOWN IN TIME ALONG RECORDS OELAYED FOR THIS VELOCITY 


DCO 305 [1TOR=20sITSPANs2 
TOR=12.5*FLOAT(ITOR) 
NEND=12 

IFLAG=0 


304 TCC=0.0 
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303 


302 


130 


IT=IT¢1 
IGR1i=I10R 


K=1 


CROSS MULTIPLY 1 SEC OF DATA FOR ALL COMBINATIONS AND 


SUM THESE AT THIS TIME 


WE Ae Be Be Be ee ee ae eK a eI 2 ee ee EK ae EA a Se a Ie ae 3K 2c 2 a a ae a a ie ak a ae a ae a 2 a Se aK aK 


DO 302 I=1+NUM 

TTCC=0.0 

ANORM1=0.0 

ANORM2=0 .0 

00 303 JJ=1+NEND 

J=JJ-1 

TT CC=TTCCHV(L( 1) s IOR1I4+J) ¥VCM(I)s TORIFJN—TO(L CI) oM(1))) 
ANORM1=ANORMI4V(L( 1) sIORI4+J) *¥*2 
ANORM2=ANORM24V (MCI) s IORI4J-TOCL (I) oM(1))) **2 

CC(LCI) oM( 1) JP=TTCC/( ANORM1 *ANORM2 ) #*0.5 


TCC=TCCHCC(L (I) oM(1I)) 


ENTER THIS VALUE INTO MATRIX 


CCORC ITs I1Ve1Z) =TCC#¥2./( AN*(AN-1.)) 
IF(L(I4+1) .EQ. LCI)) GO TO 302 
IGR1=IOR-TC(L(K) »M(K) ) 

K=K+1 


CONTINUE 
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IF CIFLAG «NE. 0) GO TO 305 
IOR=12.5*FLOAT( ITOR+1) 

IFL AG=1 

NEND=13 

GO TO 304 


305 CONTINUE 
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306 CONTINUE 
IT=1 
WRITE(6s109) A 
109 FORMAT(SXs *AZIMUTH= *,F10.22e" FOLLOWED BY V-T X-CORR 
1 PROFILE */) 
WRITE(62110) XVEE 
110 FORMAT(11Xs29F4.2/) 


ITSPAN=!I TSPAN-—20 


WRITE MATRICES 


DO 980 IT=1,ITSPAN 

WRITE(65981) ITe(CCORCIT sIVsIZ) »IV=1+IVLIM) 
981 FORMAT( 1X15 95Xs40F4.61 ) 
980 CONTINUE 

I TSPAN=ITSPAN+20 

WRITE (62983) 
983 FORMAT(////) 


107 CONTINUE 
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ITSPAN=ITSPAN-20 


STEP DOWN MATRICES IN TIME AND PICK MAXIMA ABOVE LIMIT 


ITMAX=0 

J=ITMAX4#2 

JSJ=JS+4 

IF (J «GT. ITSPAN-S) GO TO 1 

CMAX=CCOR( 10191) 

DO 91 I2Z=1.21 

00°91 1 T=J sJJ 

DO 91 IV=1eIVLIM 

IF (CCOR( ITs IVsIZ) «GE. CMAX) CMAX=CCOR(ITsIVsIZ) 
IF (CMAX .GE. CORLIM) GO TO 93 

WRITE (6952) CMAX 

FORMAT(SXs *CORREL LESS THAN LIMITs CMAX= %sF10.3) 
ITMAX=I TMAX#2 

GO TO 691 

DO 2 IZ=1.21 

DO 2 IT=JeJJ 

DO 2 IV=1eIVLIM 

IF(CCORCITsIVetIZ) «EQ. CMAX) GO TO 3 

CONTINUE 

CONTINUE 

IFC IV «LE. IVLIM-1 .AND. IV .GE. 2) GO TO 694 
WRITE(6s701) IVsIZsIT 

FORMAT (5X%e "MAX. CORREL AT MATRIX EDGEsIVeIZeIT= *% 315) 


ITMAX=I TMAX+#+2 
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GO TO 691 

IFAT AIG CEO 20 RAND. IZ .GE. 2) GO TO 693 
WRITE(6e701) IVsIZsIT 

ITMAX=ITMAX+2 

GO TO 691 


ITMAX=IT 


FIT PARABOLAS IN TIMEs VELOCITY» AND AZIMUTH 


TWOA=CCOR( IT sIV415 IZ) #CCOR( ITs I1V—-151Z)-2.*CCOR( ITs IVs IZ) 
B=(CCOR(IT sIV#191Z)—CCOR( ITs IV-191Z) )/2. — TWOAKFLOAT(IV) 
FP=-B/TWOA 

IFP=FP 


APPVEE=XVEE(IFP)+(FP-FLOAT( IFP)) *(XVEECIFRP4+1)-XVEECIFP) ) 


TWOA=CCOR(IT4191VsIZ)+CCOR( IT-15 IVeIZ)-2.*CCOR( ITs IVeIZ) 
B= (CCORCIT415IVs1IZ)—CCOR( IT—-1sIV0IZ))/2. -— TWOA*¥FLOAT(CIT) 


TP=-B/TWOA 


TWOA=CCOR( ITsIVs IZ4+1)+CCOR( ITs IVeIZ—-1)-2.*CCOR(ITsIVeIZ) 
B=(CCOR( ITs IVs IZ+1 )—CCOR( ITs IVeIZ—-1))/2. — TWOA*FLOAT(IZ) 
FP=-B/TWOA 

IFP=FP 


APP AZM=AT(IFP)+(FRP-FLOATCIFP)) *CATCIFP4+1)-ATCIFP)) 


TO=TSTRT*60.4+20. 


WRITE (62558) TO 
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134 


SS8 FORMAT(SXs°AT *%sF10.5e" SECS DCWN THE RECORD!/) 
WRITE (62560) CMAX 


S60 FORMAT(SXs "MAX. CORREL.= 'sF10.3) 


WRITE TIME AZIMUTH VELOCITY 


WRITE(6e550) APPVEEsTPsAPPAZM 
550 FORMAT(SXs *APP. VEL== "0F10.50e° TIME= %sF10.59" APP. AZM= * 
1°F10.5/) 


691 GO TO 692 
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401 CALL SKIP(NFILE-NFLST » 08) 
IFCNFILE-NFULST .LT. 0) GO TO S01 
NFUST=NFILE 
GO TG 402 

501 CALL SKIP(-150,58) 

CALL SKIP(1.0+8) 


NFULST=NF ILE 


999 CALL EXIT 


ENO 
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